This lessons can be run on EarthscapeHub. Click the button below:
This lessons can be run on EarthscapeHub. Click the button below:
Landslide susceptibility is the likelihood of a landslide occurring on the basis of local terrain conditions to estimate “where” landslides are likely to occur. Here we present a modified Data Component Use Case for Landslide Susceptibility Calculation for north Tuscany to explore how slope cohesion and soil saturation impact susceptibility. This Jupyter Notebook will first demonstrate how to use several CSDMS data components to download topography and soil datasets for an area in Tuscany, Italy then calculate the landslide susceptibility given 1) soil cohesion values and 2) soil saturation scenarios. To translate this conceptual excersice to a real-life scenario, we will later focus on a heavy rainfall event in this area on September 10-12th in 2017.
Map Source: www.mapsofworld.com
Key concepts
Generate susceptibility maps for the Tuscany region in Italy
Relationship between vegetation and soil properties to slope cohesion
Explore the influence of root chosesion (as a proxy for vegetation cover) on the lanscape's vulnerability to landslides
Relationship between cohesion, saturation, and failure probability
Factor of Safety and Susceptibility
Application Question: How might different vegetation scenarios (i.e. post wildfire) impact slope stability?
Skills
Get a general overview on how data components are implemented to get topographic and environmental data and run a landslide model
Become familiar with a simple configuration that can explore how slopes behave as we change factor od safety parameters
Make changes to root cohesion values to explore how susceptibility changes under different conditions.
Gain hands-on experience with visualizing output in Python
In this lab we will explore landslide susceptibility under different vegetation covers in the Tuscany region in Italy. To do this we modified a repository created by Gan et al. 2022 that exemplifies how models can communicate with geospatial data using the CSDMS Data components
You will run this notebook on the CSDMS Jupyterhub using the landslides kernel. You can also run it using HYDROSHARE and the CUAHSI JupyterHub (You will need to have access to both resources).
Please note that as of May 14, 2023 the Pymt library is not updated and you won't be able to run the notebook in your own PC (please contact developers for more information).
More detailed instructions for setup and environment creation the are included in the jupyter notebook from the above repository.
To use the ERA5 and Topography data componentsyou will need to create API key files to download any dataset from around the world.
Follow the instructions in CDS API Key and Open Topography API Key to get your keys.
When you have your API keys, uncomment the 2 lines in the box below to input them. Once you do it one time you can comment them again so you the code doesn't ask for the keys everytime.
# from utils import install_api_key
# install_api_key()
# import packages
import os
import warnings
import math
from numpy import inf
import numpy as np
import pandas as pd
import xarray
import xesmf as xe
import rioxarray
import cftime
from datetime import datetime
from tqdm import trange
import matplotlib.pyplot as plt
from matplotlib import colors
import imageio.v2 as imageio
from IPython.display import Video
from landlab import RasterModelGrid, imshow_grid, imshowhs_grid
from pymt.models import Topography, Era5
import matplotlib.pyplot as plt
from landlab.components import BedrockLandslider, PriorityFloodFlowRouter
from landlab.io import read_esri_ascii
from utils import regrid_data, cal_subsurface_flow_depth, cal_safety_factor
warnings.simplefilter(action='ignore', category=FutureWarning)
plt.rcParams.update({'font.size': 14})
# create folders
study_area = 'italy'
config_dir = os.path.join(os.getcwd(), 'config_files_{}'.format(study_area))
results_dir = os.path.join(os.getcwd(), 'results_{}'.format(study_area))
cache_dir = os.path.join(os.getcwd(),'cache_{}'.format(study_area))
for folder in [config_dir, results_dir, cache_dir]:
if not os.path.isdir(folder):
os.mkdir(folder)
print(folder)
Data components are used to download and visualize datasets of different nature:
Additionally:
The figure below shows the bounding box of the study area. To change the parameter settings go to 'dem_config.yaml' file.
We quickly realized that this entire area is way to big to run this code, so we will focus on in Lucca and Pistoia
Image Source: Rosi, A., Tofani, V., Tanteri, L. et al. The new landslide inventory of Tuscany (Italy) updated with PS-InSAR: geomorphological features and landslide distribution. Landslides 15, 5–19 (2018). https://doi.org/10.1007/s10346-017-0861-4
Run the following cells to:
# initialize Topography data component
dem = Topography()
dem.initialize(os.path.join(config_dir, 'dem_config.yaml'))
# get DEM variable info
var_name = dem.output_var_names[0]
var_unit = dem.var_units(var_name)
var_location = dem.var_location(var_name)
var_type = dem.var_type(var_name)
var_grid = dem.var_grid(var_name)
var_itemsize = dem.var_itemsize(var_name)
var_nbytes = dem.var_nbytes(var_name)
print('variable_name: {} \nvar_unit: {} \nvar_location: {} \nvar_type: {} \nvar_grid: {} \nvar_itemsize: {}'
'\nvar_nbytes: {} \n'. format(var_name, var_unit, var_location, var_type, var_grid, var_itemsize, var_nbytes))
variable_name: land_surface__elevation var_unit: degrees var_location: face var_type: int16 var_grid: 0 var_itemsize: 2 var_nbytes: 460800
# get DEM grid info
dem_grid_ndim = dem.grid_ndim(var_grid)
dem_grid_type = dem.grid_type(var_grid)
dem_grid_shape = dem.grid_shape(var_grid)
dem_grid_spacing = dem.grid_spacing(var_grid)
dem_grid_origin = dem.grid_origin(var_grid)
print('grid_ndim: {} \ngrid_type: {} \ngrid_shape: {} \ngrid_spacing: {} \ngrid_origin: {}'.format(
dem_grid_ndim, dem_grid_type, dem_grid_shape, dem_grid_spacing, dem_grid_origin))
grid_ndim: 2 grid_type: uniform_rectilinear grid_shape: [480 480] grid_spacing: [ 0.00083333 0.00083333] grid_origin: [ 43.85083333 10.5 ]
# get DEM variable data
dem_data = dem.get_value(var_name)
dem_data_2D = dem_data.reshape(dem_grid_shape)
# get X, Y extent for plot
min_y, min_x = dem_grid_origin
max_y = min_y + dem_grid_spacing[0]*(dem_grid_shape[0]-1)
max_x = min_x + dem_grid_spacing[1]*(dem_grid_shape[1]-1)
dy = dem_grid_spacing[0]/2
dx = dem_grid_spacing[1]/2
dem_extent = [min_x - dx, max_x + dx, min_y - dy, max_y + dy]
# plot DEM data
fig, ax = plt.subplots(1,1,figsize=(10,5))
im = ax.imshow(dem_data_2D, extent=dem_extent)
ax.title.set_text('Topography Data')
ax.set_xlabel('longitude [degree_east]')
ax.set_ylabel('latitude [degree_north]')
fig.colorbar(im,label='elevation(m)')
<matplotlib.colorbar.Colorbar at 0x7ff43cc9abb0>
# Creating RasterModelGrid of elevation data
elev = dem_data_2D.astype('float')
new_grid = RasterModelGrid(elev.shape, 90.0)
new_grid.add_field("node", "topographic__elevation", np.flipud(elev), clobber=True)
new_grid.imshow("topographic__elevation", grid_units=("deg", "deg"),
colorbar_label="Elevation (m)")
# download soil depth data
soil_raster = rioxarray.open_rasterio("https://files.isric.org/soilgrids/former/2017-03-10/data/BDRICM_M_250m_ll.tif")
soil_depth_data = soil_raster.rio.clip_box(
minx=10.50,
miny=43.50,
maxx=11.25,
maxy=44.25,
)
## plot soil depth data
soil_depth_data.plot(figsize=(6,5),cbar_kwargs={'label': 'depth(m)'})
soil_depth_data.rio.to_raster(os.path.join(cache_dir,'soil_depth.tif'))
plt.title('Soil Depth')
Text(0.5, 1.0, 'Soil Depth')
Topography data and the RasterModelGrid from Landlab to calculate the slope angle
# calculate slope angle using Topography data
model_grid = RasterModelGrid(dem_data_2D.shape,xy_spacing=(90,90))
slope = model_grid.calc_slope_at_node(elevs=dem_data) # slope in radians, 1D array
slope_angle = slope.reshape(dem_data_2D.shape) # reshape as 2D array
# # plot slope angle
fig, ax = plt.subplots(figsize=(10,5))
im=ax.imshow(slope_angle, extent=dem_extent)
cbar = fig.colorbar(im, label= 'radians')
ax.title.set_text('Slope Angle')
ax.set_xlabel('longitude [degree_east]')
ax.set_ylabel('latitude [degree_north]')
Text(0, 0.5, 'latitude [degree_north]')
All datasets downloaded or calculated above are regrid in the same grid resolution. Again for specific details refer to original notebook
# define destination grid coordinate using Topography data
dem_y = np.flip(np.arange(dem_grid_shape[0])*dem_grid_spacing[0] + dem_grid_origin[0])
dem_x = np.arange(dem_grid_shape[1])*dem_grid_spacing[1] + dem_grid_origin[1]
dest_coor = {'lon': dem_x,
'lat': dem_y}
# regrid soil depth data
soil_depth_coor = {'lon': soil_depth_data.x.values,
'lat': soil_depth_data.y.values}
soil_depth = regrid_data(soil_depth_data.values[0], soil_depth_coor, dest_coor, 'nearest_s2d')
soil_depth = soil_depth/100 # units conversion as m
# plot regridded soil depth data
fig, ax = plt.subplots(1,1,figsize=(10,5))
im = ax.imshow(soil_depth, extent=dem_extent)
ax.title.set_text('Regridded Soil Depth Data')
ax.set_xlabel('longitude [degree_east]')
ax.set_ylabel('latitude [degree_north]')
cbar = plt.colorbar(im, label='depth(m)')
Factor of safety
In geological engineering, it is common to take the ratio of the resisting stresses to driving stresses. This ratio is called the factor of safety (FS). When FS is larger than 1, the slope should be stable, while if it is below 1, the driving stress exceeds the resistance and the slope is likely to fail. FS can be calculated with the following function, and cal_safety_factor( ) is implemented based on this function.
$$ FS = \frac{(C_r + C_s)/h_s\rho_sg}{\sin\theta} + \frac{\cos\theta \tan\phi (1-\frac{h_w}{h_s}\rho_w / \rho_s)}{\sin\theta} $$where,
Susceptibility
Susceptibility is the inverse of FS. When the susceptibility is larger than 1, it means that the slope of the area is not stable and susceptible to landslide.
$$ \text{Susceptibility} = \frac{1}{FS} $$Plot Results
There are four subplots created for a list of different susceptibility values for a specified cohesioin value.
# Function for plotting raster data grid
def plot_ls(new_grid, i, cr):
# plt.figure(figsize=(6,5))
imshowhs_grid(
new_grid,
"topographic__elevation",
drape1 = np.sqrt(new_grid.at_node["susceptibility"]),
plot_type = "Drape1",
azdeg = 200,
altdeg = 20,
alpha = 0.85,
# var_name = "LS",
# var_units = r"m",
grid_units = ("(m)", "(m)"),
default_fontsize = 10,
cbar_tick_size = 10,
cmap = "hot_r",
limits = [0.5, 1.8],
ticks_km = False,
# cbar_width = "100%",
# cbar_or = "vertical",
colorbar_label_y = -55,
# add_label_bbox = True,
# bbox_to_anchor = [1.03, 0.3, 0.075, 14],
thres_drape1 = 0.01)
plt.savefig(results_dir+f'/root_cohesion_plot_{i}.png')
plt.show()
# Specify minimum and maximum root cohesion parameters for the area
min_Cr = 3000
max_Cr = 12000
# Do a list with 4 equally spaced values within that range
step_Cr = int(((max_Cr + min_Cr) - min_Cr) / 4)
root_c = [i for i in range(min_Cr, max_Cr + 1, step_Cr)] #list of values for root cohesion repoorted for area
# define mask for non-data area (water bodies, etc)
mask = (slope_angle == 0) & (soil_depth > 2.0)
# we assume constant saturation so that relative wetness is one
subsurface_flow_depth = soil_depth
# plot susceptibility
# fig = plt.figure(figsize=(16,10))
# fig.suptitle("Susceptibility of landslide for different root cohesion values")
#change plot grid size if you want to explore more parameters in the above list
nrows, ncols = 2,2
# calculate FS for each root cohesion value
for i,cr in enumerate(root_c):
safety_factor = cal_safety_factor(slope_angle, subsurface_flow_depth, soil_depth,
root_cohesion=cr, soil_cohesion=5000, soil_bulk_density=1300,
soil_internal_friction_angle=35)
# calculate susceptibility
susceptibility = 1.0 / safety_factor
susceptibility = np.where(~mask, susceptibility, np.nan) #water bodies
susceptibility = np.where(susceptibility > 0.5, susceptibility, 0)
susp = new_grid.add_field("node", "susceptibility", np.flipud(susceptibility), clobber=True)
ax = fig.add_subplot(nrows, ncols, i+1)
plt.title('Cr = '+str(cr)+" (Pa kg/ms$^{2}$)", loc = 'center', fontdict={'fontsize': 12})
plot_ls(new_grid, i, cr)
# im_data = susceptibility
###### Comment the 2 above lines and uncomment below to check the differences""
# if cr == 5000:
# ax = fig.add_subplot(nrows, ncols, 1)
# im_data = susceptibility
# reference_sus = np.copy(susceptibility)
# else:
# ax = fig.add_subplot(nrows, ncols, i+1)
# im_data = susceptibility - reference_sus
# im_sus = ax.imshow(im_data, cmap='BrBG_r', extent=dem_extent)
# plt.colorbar(im_sus, ax=ax)
# im_sus.set_clim(0,1.1)
# plt.tight_layout()
Remember that a susceptibility close to 1 means the landscape is prone to landslides and when it's close to zero it is stable.
Now, we will see how our maps change as we vary two parameters in the Factory of Safety equation, root cohesion and soil saturation
Here we look at three (3) different saturation level scenarios (0%, 50%, and 100%)
Again, remember that a susceptibility close to 1 means the landscape is prone to landslides and when it's close to zero it is stable.
# Specify minimum and maximum root cohesion parameters for the area
min_Cr=3000
max_Cr=12000
# Do a list with 4 equally spaced values within that range
step_Cr=int(((max_Cr+min_Cr)-min_Cr)/3)
root_c=[i for i in range(min_Cr,max_Cr+1,step_Cr)] #list of values for root cohesion repoorted for area
#root_c=[3000, 8000, 12000]
#sats=[0.1, 0.5, 1.0]
# define mask for non-data area (water bodies, etc)
mask = (slope_angle==0)&(soil_depth>2.0)
# plot susceptibility
fig = plt.figure(figsize=(12,10))
################# Want to put x axis label at top, want tick marks centered on the columns #####################
################# Want y axis ticks centered on rows ###########################################################
# fig.suptitle("Susceptibility of landslide", y=0.95) #for different soil saturations and root cohesion values")
#change plot grid size if you want to explore more parameters in the above list
nrows, ncols = 3,3 #rows are cohesion values, columns are saturation
gs = fig.add_gridspec(nrows=nrows+2, ncols=ncols+1, width_ratios=[0.1, 1, 1, 1], height_ratios=[0.1, 1, 1, 1, 0.1], wspace=0.2,hspace=0.3)
xlab_ax = fig.add_subplot(gs[0,:])
ylab_ax = fig.add_subplot(gs[1:-1,0])
xlab_ax.text(x=0.5, y=0.05, ha='center', s='Saturation %')
ylab_ax.text(x=0.05, y=0.5, rotation=90, va='center', s='Root Cohesion (Pa kg/ms^2)')
xlab_ax.axis('off')
ylab_ax.axis('off')
w=0
percent = np.arange(0,11,5)
sat_depth = [p/10 * soil_depth for p in percent]
for i, cr in enumerate(root_c, start=1):
for j, subsurface_flow_depth in enumerate(sat_depth, start=1): # Only interested in 0%, 50%, and 100% Saturation
# if j==0:
# subsurface_flow_depth = 0
safety_factor = cal_safety_factor(slope_angle, subsurface_flow_depth, soil_depth,
root_cohesion=cr, soil_cohesion=5000, soil_bulk_density=1300,
soil_internal_friction_angle=35)
# calculate susceptibility
susceptibility = 1.0 / safety_factor
susceptibility = np.where(~mask, susceptibility, np.nan) #water bodies
# Changing the index
w=w+1
# plot
# ax = fig.add_subplot(nrows,ncols,w)
ax = fig.add_subplot(gs[i,j], box_aspect=1)
im_data = susceptibility
ax.tick_params(axis='both', which='major', labelsize=8)
ax.tick_params(axis='both', which='minor', labelsize=8)
im_sus = ax.imshow(im_data, cmap='BrBG_r', extent=dem_extent)
if j == 1:
ax.set_ylabel(cr)
if i == 1:
ax.set_xlabel(10*percent[j-1])
ax.xaxis.set_label_position('top')
# plt.colorbar(im_sus, ax=ax)
im_sus.set_clim(0,1.0)
# plt.close(fig)
cb_ax = fig.add_subplot(gs[-1, :])
cbar = fig.colorbar(im_sus, cax=cb_ax, location='bottom', orientation='horizontal')
cbar.set_label('Susceptibility of Landslide')
# im_sus.set_clim(0,1.1)
# fig.savefig(os.path.join(results_dir, 'sus_root_{}.png'.format(j)))
This process is similiar to what we did for the original root cohesion comparison
# Function for plotting raster data grid
def plot_ls(new_grid, i, j):
# plt.figure(figsize=(6,5))
imshowhs_grid(
new_grid,
"topographic__elevation",
drape1 = np.sqrt(new_grid.at_node["susceptibility"]),
plot_type = "Drape1",
azdeg = 200,
altdeg = 20,
alpha = 0.85,
# var_name = "LS",
# var_units = r"m",
grid_units = ("(m)", "(m)"),
default_fontsize = 10,
cbar_tick_size = 10,
cmap = "hot_r",
limits = [0.5, 1.8],
ticks_km = False,
# cbar_width = "100%",
# cbar_or = "vertical",
colorbar_label_y = -55,
# add_label_bbox = True,
# bbox_to_anchor = [1.03, 0.3, 0.075, 14],
thres_drape1 = 0.01)
plt.savefig(results_dir+f'/sus_root_plot_{i}_{j}.png')
plt.show()
# Specify minimum and maximum root cohesion parameters for the area
min_Cr=3000
max_Cr=12000
# Do a list with 4 equally spaced values within that range
step_Cr=int(((max_Cr+min_Cr)-min_Cr)/3)
root_c=[i for i in range(min_Cr,max_Cr+1,step_Cr)] #list of values for root cohesion repoorted for area
#root_c=[3000, 8000, 12000]
#sats=[0.1, 0.5, 1.0]
# define mask for non-data area (water bodies, etc)
mask = (slope_angle==0)&(soil_depth>2.0)
# plot susceptibility
# fig = plt.figure(figsize=(16,16))
# plt.xlabel('Saturation %')
# plt.ylabel('Root Cohesion (Pa kg/ms^2)')
# cb_ax = fig.add_axes([1.02, 0.1, 0.02, 0.8])
# cbar = fig.colorbar(im_sus, cax=cb_ax)
# im_sus.set_clim(0,1.1)
################# Want to put x axis label at top, want tick marks centered on the columns #####################
################# Want y axis ticks centered on rows ###########################################################
fig.suptitle("Susceptibility of landslide") #for different soil saturations and root cohesion values")
#change plot grid size if you want to explore more parameters in the above list
nrows, ncols = 3,3 #rows are cohesion values, columns are saturation
w=0
for i,cr in enumerate(root_c):
for j in range(0,11,5): # Only interested in 0%, 50%, and 100% Saturation
if j==0:
subsurface_flow_depth = 0
subsurface_flow_depth = soil_depth*j/10
safety_factor = cal_safety_factor(slope_angle, subsurface_flow_depth, soil_depth,
root_cohesion=cr, soil_cohesion=5000, soil_bulk_density=1300,
soil_internal_friction_angle=35)
# calculate susceptibility
susceptibility = 1.0 / safety_factor
susceptibility = np.where(~mask, susceptibility, np.nan) #water bodies
susceptibility = np.where(susceptibility > 0.5, susceptibility, 0)
susp = new_grid.add_field("node", "susceptibility", np.flipud(susceptibility), clobber=True)
# ax = fig.add_subplot(nrows, ncols, i+1)
plt.title('Cr = '+str(cr)+" (Pa kg/ms$^{2}$)", loc = 'center', fontdict={'fontsize': 12})
plot_ls(new_grid, i, j)
If the model is run in its complete form:
volumetric soil water data will be used for calculating the susceptibility
precipitation data will help visualize
The 'era5_config.yaml' file includes the parameter settings of this data component.
Run the following cells to:
# initialize ERA5 data component
era5 = Era5()
era5.initialize(os.path.join(config_dir,'era5_config.yaml'))
2023-05-17 04:42:09,421 INFO Welcome to the CDS 2023-05-17 04:42:09,421 INFO Sending request to https://cds.climate.copernicus.eu/api/v2/resources/reanalysis-era5-single-levels 2023-05-17 04:42:09,809 INFO Request is queued 2023-05-17 04:42:10,982 INFO Request is running 2023-05-17 04:42:31,644 INFO Request is completed 2023-05-17 04:42:31,645 INFO Downloading https://download-0007-clone.copernicus-climate.eu/cache-compute-0007/cache/data8/adaptor.mars.internal-1684298544.1200817-10618-3-07c501a3-bafd-41b3-87de-3675db8f0d75.nc to ./cache_italy/era5.nc (6.2K) 2023-05-17 04:42:32,488 INFO Download rate 7.4K/s
# get ERA5 variable info
for var_name in era5.output_var_names:
var_unit = era5.var_units(var_name)
var_location = era5.var_location(var_name)
var_type = era5.var_type(var_name)
var_grid = era5.var_grid(var_name)
var_itemsize = era5.var_itemsize(var_name)
var_nbytes = era5.var_nbytes(var_name)
# print('variable_name: {} \nvar_unit: {} \nvar_location: {} \nvar_type: {} \nvar_grid: {} \nvar_itemsize: {}'
# '\nvar_nbytes: {} \n'. format(var_name, var_unit, var_location, var_type, var_grid, var_itemsize, var_nbytes))
# get ERA5 grid info
era5_grid_ndim = era5.grid_ndim(var_grid)
era5_grid_type = era5.grid_type(var_grid)
era5_grid_shape = era5.grid_shape(var_grid)
era5_grid_spacing = era5.grid_spacing(var_grid)
era5_grid_origin = era5.grid_origin(var_grid)
# print('grid_ndim: {} \ngrid_type: {} \ngrid_shape: {} \ngrid_spacing: {} \ngrid_origin: {}'.format(
# era5_grid_ndim, era5_grid_type, era5_grid_shape, era5_grid_spacing, era5_grid_origin))
# get ERA5 time info
era5_start_time = era5.start_time
era5_end_time = era5.end_time
era5_time_step = era5.time_step
era5_time_unit = era5.time_units
era5_time_steps = int((era5_end_time - era5_start_time)/era5_time_step) + 1
# print('start_time:{} \nend_time:{} \ntime_step:{} \ntime_unit:{} \ntime_steps:{}'.format(
# era5_start_time, era5_end_time, era5_time_step, era5_time_unit, era5_time_steps))
# get ERA5 variables data and plot (at the first time step)
fig = plt.figure(figsize=(18,16))
nrows, ncols = 3, 2
i = 1
for var_name in era5.output_var_names:
ax = fig.add_subplot(nrows, ncols, i)
var_unit = era5.var_units(var_name)
# get variable data
era5_data = era5.get_value(var_name)
era5_data_2D = era5_data.reshape(era5_grid_shape)
# get X, Y extent for plot
min_y, min_x = era5_grid_origin
max_y = min_y + era5_grid_spacing[0]*(era5_grid_shape[0]-1)
max_x = min_x + era5_grid_spacing[1]*(era5_grid_shape[1]-1)
dy = era5_grid_spacing[0]/2
dx = era5_grid_spacing[1]/2
era5_extent = [min_x - dx, max_x + dx, min_y - dy, max_y + dy]
# plot data
im = ax.imshow(era5_data_2D, extent=era5_extent, cmap='Blues')
ax.title.set_text('{} ({})'.format(var_name,var_unit ))
ax.set_xlabel('longitude [degree_east]')
ax.set_ylabel('latitude [degree_north]')
cbar = plt.colorbar(im, ax=ax)
plt.tight_layout()
i += 1
